RaceID progenitor

library(RaceID)
library(Seurat)

library(dplyr)

library(ggplot2)
library(patchwork)

Attaching SeuratObject


Attaching package: ‘dplyr’


The following objects are masked from ‘package:stats’:

    filter, lag


The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

# Set working directory to project root
setwd("/research/peer/fdeckert/FD20200109SPLENO")

# Source files
source("plotting_global.R")

Loading required package: RColorBrewer

Parameter settings

# Seurat files 
so_file <- "data/object/seurat.rds"

# Plotting Theme
ggplot2::theme_set(theme_global_set()) # From project global source()

Import Seurat

so <- readRDS(so_file)

so <- subset(so, subset = nFeature_RNA <= 3000 & nCount_RNA <= 15000)

Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from rna_umap_ to rnaumap_”

Split data bei treatment and tissue group

so_1 <- subset(so, subset = treatment == "NaCl" & tissue == "Progenitor" & rna_snn_res.0.8 != 18)
so_2 <- subset(so, subset = treatment == "CpG" & tissue == "Progenitor" & rna_snn_res.0.8 != 18)

Get differently expressed genes between condition

# library(Seurat)

# Idents(so_1) <- "sample_name"
# deg_1 <- FindMarkers(so_1, ident.1 = "Progenitor_NaCl_Rep1", ident.2 = "Progenitor_NaCl_Rep2")

# Idents(so_2) <- "sample_name"
# deg_2 <- FindMarkers(so_2, ident.1 = "Progenitor_CpG_Rep1", ident.2 = "Progenitor_CpG_Rep2")

Create SCseq Object

sc_1 <- SCseq(expdata = GetAssayData(so_1, slot = "counts", assay = "RNA"))
sc_2 <- SCseq(expdata = GetAssayData(so_2, slot = "counts", assay = "RNA"))

SCseq filter data

Filterdata can be run with LBatch and bmode to remove genes associated with batch effects. VarID workflow also allows regression of batch effect in pruneKnn with regNB. in that case the regNB parameter of the compNoise function needs to be TRUE as well. Additional variables can be regressed out using the regVar variable with a data.frame. The rownames correspond to cell id and columns to values e.g. cell cycle.

# Filter
sc_1 <- filterdata(sc_1, mintotal = 2500, minexpr = 5, minnumber = 5, CGenes = rownames(so_1)[grep("^(mt)", rownames(so_1))], FGenes = NULL)
sc_2 <- filterdata(sc_2, mintotal = 2500, minexpr = 5, minnumber = 5, CGenes = rownames(so_2)[grep("^(mt)", rownames(so_2))], FGenes = NULL)

Filtered expression matrix

expdata_1 <- getExpData(sc_1)
expdata_2 <- getExpData(sc_2)

Create variables for regression

# Batch 
batch_1 <- so_1$sample_name
names(batch_1) <- colnames(so_1)
batch_1 <- batch_1[names(batch_1) %in% colnames(expdata_1)]

batch_2 <- so_2$sample_name
names(batch_2) <- colnames(so_2)
batch_2 <- batch_2[names(batch_2) %in% colnames(expdata_2)]

Prune KNN and compute cluster with transition probabilities

res_1   <- pruneKnn(expdata_1, large = TRUE, regNB = TRUE, batch = batch_1, knn = 10, alpha = NULL, no_cores = NULL, seed = 12345)
cl_1    <- graphCluster(res_1, pvalue = 0.01)
probs_1 <- transitionProbs(res_1, cl_1)

res_2   <- pruneKnn(expdata_2, large = TRUE, regNB = TRUE, batch = batch_2, knn = 10, alpha = NULL, no_cores = NULL, seed = 12345)
cl_2    <- graphCluster(res_2, pvalue = 0.01)
probs_2 <- transitionProbs(res_2, cl_2)

Compute noise from corrected variance

noise_1 <- compNoise(expdata_1, res_1, regNB = TRUE, batch = batch_1, pvalue = 0.01, no_cores = NULL)
noise_2 <- compNoise(expdata_2, res_2, regNB = TRUE, batch = batch_2, pvalue = 0.01, no_cores = NULL)

Update SCseq object

sc_1 <- updateSC(sc_1, res = res_1, cl = cl_1, noise = noise_1, flo = .1)
sc_2 <- updateSC(sc_2, res = res_2, cl = cl_2, noise = noise_2, flo = .1)

Save output

saveRDS(sc_1, "data/object/components/raceid_progenitor_mt_nacl.rds")
saveRDS(sc_2, "data/object/components/raceid_progenitor_mt_cpg.rds")

saveRDS(res_1, "data/object/components/raceid_progenitor_mt_nacl_res.rds")
saveRDS(res_2, "data/object/components/raceid_progenitor_mt_cpg_res.rds")

saveRDS(cl_1, "data/object/components/raceid_progenitor_mt_nacl_cl.rds")
saveRDS(cl_2, "data/object/components/raceid_progenitor_mt_cpg_cl.rds")

saveRDS(probs_1, "data/object/components/raceid_progenitor_mt_nacl_probs.rds")
saveRDS(probs_2, "data/object/components/raceid_progenitor_mt_cpg_probs.rds")

# sc_1 <- readRDS("data/object/components/raceid_progenitor_mt_nacl.rds")
# sc_2 <- readRDS("data/object/components/raceid_progenitor_mt_cpg.rds")

# res_1 <- readRDS("data/object/components/raceid_progenitor_mt_nacl_res.rds")
# res_2 <- readRDS("data/object/components/raceid_progenitor_mt_cpg_res.rds")

# cl_1 <- readRDS("data/object/components/raceid_progenitor_mt_nacl_cl.rds")
# cl_2 <- readRDS("data/object/components/raceid_progenitor_mt_cpg_cl.rds")

# probs_1 <- readRDS("data/object/components/raceid_progenitor_mt_nacl_probs.rds")
# probs_2 <- readRDS("data/object/components/raceid_progenitor_mt_cpg_probs.rds")

Dimensional reduction

sc_1 <- compumap(sc_1, spread = 3)
sc_2 <- compumap(sc_2, spread = 3)

Plot

options(repr.plot.width = 5, repr.plot.height = 5)
umap_1 <- plotmap(sc_1, um = TRUE)
umap_2 <- plotmap(sc_2, um = TRUE)

## map of transition probabilities
options(repr.plot.width = 5, repr.plot.height = 5)
plotTrProbs(sc_1, probs_1, tp = .5, prthr = 0.01, cthr = 0, um = TRUE)
plotTrProbs(sc_2, probs_2, tp = .5, prthr = 0.01, cthr = 0, um = TRUE)

Combine VarID results Seurat object

so_1 <- subset(so_1, cells = rownames(sc_1@umap))
so_2 <- subset(so_2, cells = rownames(sc_2@umap))

all(colnames(so_1) == rownames(sc_1@umap))
all(colnames(so_2) == rownames(sc_2@umap))

TRUE

TRUE

so_1$cluster <- sc_1@cluster$kpart
so_2$cluster <- sc_2@cluster$kpart

so_1@reductions$varid_umap <- CreateDimReducObject(embeddings = as.matrix(sc_1@umap), assay = "RNA", key = "varid_umap")
so_2@reductions$varid_umap <- CreateDimReducObject(embeddings = as.matrix(sc_2@umap), assay = "RNA", key = "varid_umap")

Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from varid_umap to varidumap_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to varidumap_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from varid_umap to varidumap_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to varidumap_”

UMAP reduction plots

reduction <- "varid_umap"

options(repr.plot.width = 20, repr.plot.height = 10)

dplot_1 <- DimPlot(so_1, reduction = reduction, group.by = "sample_name", label = FALSE) & 
    ggtitle("NaCl") & xlab("UMAP 1") & ylab("UMAP 2") &
    theme(aspect.ratio = 1, legend.position = "bottom") & 
    guides(color = guide_legend(ncol = 1, override.aes = list(size = 2)))

dplot_2 <- DimPlot(so_1, reduction = reduction, group.by = "cluster", label = TRUE) & 
    ggtitle("NaCl") & xlab("UMAP 1") & ylab("UMAP 2") &
    theme(aspect.ratio = 1, legend.position = "none") & 
    guides(color = guide_legend(ncol = 3, override.aes = list(size = 2)))

dplot_3 <- DimPlot(so_1, reduction = reduction, group.by = "main_labels", label = FALSE) & 
    ggtitle("NaCl") & xlab("UMAP 1") & ylab("UMAP 2") &
    theme(aspect.ratio = 1, legend.position = "bottom") & 
    scale_color_manual(values = color$main_labels_haemosphere, breaks = names(color$main_labels_haemosphere), na.value = "dark gray") & 
    guides(color = guide_legend(ncol = 3, override.aes = list(size = 2)))

dplot_4 <- DimPlot(so_1, reduction = reduction, group.by = "fine_labels", label = FALSE) & 
    ggtitle("NaCl") & xlab("UMAP 1") & ylab("UMAP 2") &
    theme(aspect.ratio = 1, legend.position = "bottom") & 
    scale_color_manual(values = color$fine_labels_haemosphere, breaks = names(color$fine_labels_haemosphere), na.value = "dark gray") & 
    guides(color = guide_legend(ncol = 3, override.aes = list(size = 2)))

dplot_5 <- DimPlot(so_2, reduction = reduction, group.by = "sample_name", label = FALSE) & 
    ggtitle("CpG") & xlab("UMAP 1") & ylab("UMAP 2") &
    theme(aspect.ratio = 1, legend.position = "bottom") & 
    guides(color = guide_legend(ncol = 1, override.aes = list(size = 2)))

dplot_6 <- DimPlot(so_2, reduction = reduction, group.by = "cluster", label = TRUE) & 
    ggtitle("CpG") & xlab("UMAP 1") & ylab("UMAP 2") &
    theme(aspect.ratio = 1, legend.position = "none") & 
    guides(color = guide_legend(ncol = 3, override.aes = list(size = 2)))

dplot_7 <- DimPlot(so_2, reduction = reduction, group.by = "main_labels", label = FALSE) & 
    ggtitle("CpG") & xlab("UMAP 1") & ylab("UMAP 2") &
    theme(aspect.ratio = 1, legend.position = "bottom") & 
    scale_color_manual(values = color$main_labels_haemosphere, breaks = names(color$main_labels_haemosphere), na.value = "dark gray") & 
    guides(color = guide_legend(ncol = 3, override.aes = list(size = 2)))

dplot_8 <- DimPlot(so_2, reduction = reduction, group.by = "fine_labels", label = FALSE) & 
    ggtitle("CpG") & xlab("UMAP 1") & ylab("UMAP 2") &
    theme(aspect.ratio = 1, legend.position = "bottom") & 
    scale_color_manual(values = color$fine_labels_haemosphere, breaks = names(color$fine_labels_haemosphere), na.value = "dark gray") & 
    guides(color = guide_legend(ncol = 3, override.aes = list(size = 2)))


dplot <- dplot_1 + dplot_2 + dplot_3 + dplot_4 + dplot_5 + dplot_6 + dplot_7 + dplot_8 + plot_layout(ncol = 4)
dplot

fplot_1 <- FeaturePlot(so_1, reduction = reduction, features = "nCount_RNA") & theme(aspect.ratio = 1)
fplot_2 <- FeaturePlot(so_1, reduction = reduction, features = "nFeature_RNA") & theme(aspect.ratio = 1)
fplot_3 <- FeaturePlot(so_1, reduction = reduction, features = "fine_delta_score") & theme(aspect.ratio = 1)

fplot_4 <- FeaturePlot(so_2, reduction = reduction, features = "nCount_RNA") & theme(aspect.ratio = 1)
fplot_5 <- FeaturePlot(so_2, reduction = reduction, features = "nFeature_RNA") & theme(aspect.ratio = 1)
fplot_6 <- FeaturePlot(so_2, reduction = reduction, features = "fine_delta_score") & theme(aspect.ratio = 1)

options(repr.plot.width = 15, repr.plot.height = 10)
fplot <- fplot_1 + fplot_2 + fplot_3 +  fplot_4 + fplot_5 + fplot_6 + plot_layout(ncol = 3)
fplot

fplot_1 <- FeaturePlot(so_1, reduction = reduction, features = "pMt_RNA") & theme(aspect.ratio = 1)
fplot_2 <- FeaturePlot(so_1, reduction = reduction, features = "pRp_RNA") & theme(aspect.ratio = 1)
fplot_3 <- FeaturePlot(so_1, reduction = reduction, features = "pHb_RNA") & theme(aspect.ratio = 1)

fplot_4 <- FeaturePlot(so_2, reduction = reduction, features = "pMt_RNA") & theme(aspect.ratio = 1)
fplot_5 <- FeaturePlot(so_2, reduction = reduction, features = "pRp_RNA") & theme(aspect.ratio = 1)
fplot_6 <- FeaturePlot(so_2, reduction = reduction, features = "pHb_RNA") & theme(aspect.ratio = 1)

options(repr.plot.width = 15, repr.plot.height = 10)
fplot <- fplot_1 + fplot_2 + fplot_3 +  fplot_4 + fplot_5 + fplot_6 + plot_layout(ncol = 3)
fplot

fplot_1 <- FeaturePlot(so_1, reduction = reduction, features = "msG0_RNA1") & theme(aspect.ratio = 1)
fplot_2 <- FeaturePlot(so_1, reduction = reduction, features = "msInt_RNA1") & theme(aspect.ratio = 1)
fplot_3 <- FeaturePlot(so_1, reduction = reduction, features = "msCC_diff_RNA") & theme(aspect.ratio = 1)

fplot_4 <- FeaturePlot(so_2, reduction = reduction, features = "msG0_RNA1") & theme(aspect.ratio = 1)
fplot_5 <- FeaturePlot(so_2, reduction = reduction, features = "msInt_RNA1") & theme(aspect.ratio = 1)
fplot_6 <- FeaturePlot(so_2, reduction = reduction, features = "msCC_diff_RNA") & theme(aspect.ratio = 1)

options(repr.plot.width = 15, repr.plot.height = 10)
fplot <- fplot_1 + fplot_2 + fplot_3 +  fplot_4 + fplot_5 + fplot_6 + plot_layout(ncol = 3)
fplot

Cluster label frequency NaCl

cluster_main_labels <- ggplot(so_1@meta.data, aes(x = factor(cluster), fill = factor(main_labels, levels = names(color$main_labels))))+ 
  geom_bar(stat = "count", position = "fill") + 
  scale_fill_manual(values = color$main_labels_haemosphere, breaks = names(color$main_labels_haemosphere), na.value = "dark gray") +
  ggtitle("Cluster frequency progenitor NaCl") + xlab("Cluster") + ylab("Cell frequency") + 
  theme(legend.position = "bottom") + 
    guides(fill = guide_legend(ncol = 5, override.aes = list(size = 2), title = "Main labels"))

all_main_labels <- ggplot(so_1@meta.data, aes(x = "NaCl", fill = factor(main_labels, levels = names(color$main_labels))))+ 
  geom_bar(stat = "count", position = "fill", width = 1) + 
  scale_fill_manual(values = color$main_labels_haemosphere, breaks = names(color$main_labels_haemosphere), na.value = "dark gray") +
  ggtitle("Sort frequency progenitor NaCl") + xlab("") + ylab("Cell frequency") + 
  theme(legend.position = "bottom") + 
    guides(fill = guide_legend(ncol = 5, override.aes = list(size = 2), title = "Main labels"))

cluster_fine_labels <- ggplot(so_1@meta.data, aes(x = factor(cluster), fill = factor(fine_labels, levels = names(color$fine_labels)))) + 
  geom_bar(stat = "count", position = "fill") + 
  scale_fill_manual(values = color$fine_labels_haemosphere, breaks = names(color$fine_labels_haemosphere), na.value = "dark gray") +
  ggtitle("Cluster frequency progenitor NaCl") + xlab("Cluster") + ylab("Cell frequency") + 
  theme(legend.position = "bottom") + 
  guides(fill = guide_legend(ncol = 5, override.aes = list(size = 2), title = "Fine labels"))

all_fine_labels <- ggplot(so_1@meta.data, aes(x = "NaCl", fill = factor(fine_labels, levels = names(color$fine_labels)))) + 
  geom_bar(stat = "count", position = "fill", width = 1) + 
  scale_fill_manual(values = color$fine_labels_haemosphere, breaks = names(color$fine_labels_haemosphere), na.value = "dark gray") +
  ggtitle("Cluster frequency progenitor NaCl") + xlab("Cluster") + ylab("Cell frequency") + 
  theme(legend.position = "bottom") + 
  guides(fill = guide_legend(ncol = 5, override.aes = list(size = 2), title = "Fine labels"))

options(repr.plot.width = 20, repr.plot.height = 5)
cluster_main_labels + all_main_labels + cluster_fine_labels + all_fine_labels + plot_layout(width = c(1, 1/11, 1, 1/11), ncol = 4, guides = 'collect') & theme(legend.position = 'bottom')

Cluster label frequency CpG

cluster_main_labels <- ggplot(so_2@meta.data, aes(x = factor(cluster), fill = factor(main_labels, levels = names(color$main_labels))))+ 
  geom_bar(stat = "count", position = "fill") + 
  scale_fill_manual(values = color$main_labels_haemosphere, breaks = names(color$main_labels_haemosphere), na.value = "dark gray") +
  ggtitle("Cluster frequency progenitor NaCl") + xlab("Cluster") + ylab("Cell frequency") + 
  theme(legend.position = "bottom") + 
    guides(fill = guide_legend(ncol = 5, override.aes = list(size = 2), title = "Main labels"))

all_main_labels <- ggplot(so_2@meta.data, aes(x = "NaCl", fill = factor(main_labels, levels = names(color$main_labels))))+ 
  geom_bar(stat = "count", position = "fill", width = 1) + 
  scale_fill_manual(values = color$main_labels_haemosphere, breaks = names(color$main_labels_haemosphere), na.value = "dark gray") +
  ggtitle("Sort frequency progenitor NaCl") + xlab("") + ylab("Cell frequency") + 
  theme(legend.position = "bottom") + 
    guides(fill = guide_legend(ncol = 5, override.aes = list(size = 2), title = "Main labels"))

cluster_fine_labels <- ggplot(so_2@meta.data, aes(x = factor(cluster), fill = factor(fine_labels, levels = names(color$fine_labels)))) + 
  geom_bar(stat = "count", position = "fill") + 
  scale_fill_manual(values = color$fine_labels_haemosphere, breaks = names(color$fine_labels_haemosphere), na.value = "dark gray") +
  ggtitle("Cluster frequency progenitor NaCl") + xlab("Cluster") + ylab("Cell frequency") + 
  theme(legend.position = "bottom") + 
  guides(fill = guide_legend(ncol = 5, override.aes = list(size = 2), title = "Fine labels"))

all_fine_labels <- ggplot(so_2@meta.data, aes(x = "NaCl", fill = factor(fine_labels, levels = names(color$fine_labels)))) + 
  geom_bar(stat = "count", position = "fill", width = 1) + 
  scale_fill_manual(values = color$fine_labels_haemosphere, breaks = names(color$fine_labels_haemosphere), na.value = "dark gray") +
  ggtitle("Cluster frequency progenitor NaCl") + xlab("Cluster") + ylab("Cell frequency") + 
  theme(legend.position = "bottom") + 
  guides(fill = guide_legend(ncol = 5, override.aes = list(size = 2), title = "Fine labels"))

options(repr.plot.width = 20, repr.plot.height = 5)
cluster_main_labels + all_main_labels + cluster_fine_labels + all_fine_labels + plot_layout(width = c(1, 1/11, 1, 1/11), ncol = 4, guides = 'collect') & theme(legend.position = 'bottom')